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C^ Abstract. Mesenchymal motion describes the movement of cells in biological 

G tissues formed by fibre networks. An important example is the migration of 

I "^ I tumour cells through collagen networks during the process of metastasis forma- 
tion. We investigate the mesenchymal motion model proposed by T. Hillen in 

/»f~\ [14] in higher dimensions. We formulate the problem as an evolution equation 

^ in a Banach space of measure-valued functions and use methods from semi- 

^Tv group theory to show the global existence of classical solutions. We investigate 

^4- steady states of the model and show that patterns of network type exist as 

<-0 steady states. For the case of constant fibre distribution, we find an explicit 

^vj solution and we prove the convergence to the parabolic limit. 

^D 1. Introduction. Friedl and collaborators [11] observed mesenchymal tumour cells 

as they move in a field of collagen fibres and change their velocities according to 
the local orientation of the fibres. At the same time, the cells also remodel the 
^ fibres, primarily through expression of matrix-degrading enzymes (proteases) that 

k>< cut selected fibres. In [14], the author introduced a mathematical model for this 

^ process of mesenchymal cell movement in fibrous tissues. Recent analysis of this 

Cu and similar models [14, 21, 4, 5, 25] revealed the existence of biologically meaningful 

measure valued solutions, which correspond to tissue and cell alignment. Hence a 
sophisticated existence theory is needed. In this paper we will formulate the mes- 
enchymal transport model proposed in [14] as a semilinear evolution equation in 
a Banach space of measure- valued functions. We apply classical theory of semi- 
groups of operators and a Banach Fixed Point argument to show well-posedness of 
the problem (Section 3.1). With the correct theoretical framework in place, we are 
then able to classify possible steady states, whereby we introfduce a new notation 
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of pointwise steady states, which are meant to resemble the network patterns which 
were observed numericaUy in [21]. Moreover, we rigorously study the parabolic 
limit (diffusion limit) of the kinetic model in the measure-valued context. We show 
convergence to the diffusion limit for constant fibre distribution in Section 5.1. 

The existence theory here employs a mild solution formulation which is based 
on a variation of constant formula. The solutions are functions in L^ in space and 
measures in velocity. It turns out that this definition is too "weak" in the sense that 
it does not provide a nice representation of the global network patterns observed 
numerically. Hence here we introduce a sub-class of steady states, which we call 
pointwise steady states. First of all, we show that pointwise steady states do exist. 
Secondly, pointwise steady states allow for a representation of network patterns. 
Our results include a classification of possible network intersections. 

In the model proposed in [14], undirected and directed tissues were distinguished. 
In undirected tissues (e.g. collagen), fibres are symmetric and both directions are 
identical, a situation that somewhat resembles a nematic liquid crystal [24]. In 
directed tissues, fibres are asymmetric and the two ends can be distinguished. From 
the mathematical point of view, which we adopt in the present paper, both cases 
are completely analogous. Hence we focus on the case of undirected tissues. We 
refer to [14] for the biological assumptions and the detailed mathematical derivation 
of the model. 

The model studied here is specifically designed for mesenchymal cell movement 
in network tissues via contact guidance and degradation of the extracellular matrix 
(ECM). Painter [21] has extended this model in various directions. His model vari- 
ations allow (i) to choose between amoeboid and mesenchymal motion, (ii) to place 
different weights between diffusive movement and movement by contact guidance, 
(iii) to include ECM degradation as well as production, (iv) to include ECM remod- 
elling or lack thereof, (v) to study focussed protease release at the cell tip versus 
unfocussed ECM degradation via a diffusible proteolytic enzyme. Many of these 
modifications lead to the same pattern formation properties as observed for the 
initial model. All of these modifications show the same mathematical challenges, 
namely the description of aligned tissue as weak solutions and orientation driven 
instabilities. Hence we believe that the results which we present here are represen- 
tative for a large class of kinetic models for cell movement in tissues and they can 
be generalized to many other cases. 

In [14], the techniques of moment closure, parabolic and hydrodynamic scaling 
were used to study the macroscopic limits of the system that we later restate in 
equation (1). The resulting macroscopic models have the form of drift-diffusion 
equations where the mean drift velocity is given by the mean orientation of the 
tissue and the diffusion tensor is given by the variance-covariance matrix of the tissue 
orientations. Model (1) has been extended in [4, 5] to include cell-cell interactions 
and chemotaxis. The corresponding diffusion limit was formally obtained in these 
papers. 

In the case of chemotaxis, a system of a transport equation for the cell motion 
coupled to a parabolic or elliptic equation for the chemical signal was studied by 
Alt [1], Chalub et al. [3] and Hwang et al. [15, 16]. Local and global existence 
of solutions were studied and the macroscopic limits were proved rigorously in [3, 
15, 16]. However, these authors assumed the existence of an equilibrium velocity 
distribution for cells that is in L°°{V) where V denotes the space of velocities. For 
the mesenchymal motion model, it is necessary to allow for complete alignments 
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of either fibres or cells, corresponding to Dirac measures on V or the space of 
directions, the unit sphere §"~^. In particular, assumption (AO) in paper [3] does 
not apply here and hence their respective results can not be applied directly to the 
mesenchymal motion model. 

In Section 2 we formulate the model and we introduce suitable function spaces 
and operators. Our first main result on global existence of measure- valued solutions 
is given in Section 3. In Section 4 we present a definition and classification of 
pointwise steady states. In Section 5 we assume that the fibre density g is a given 
function of x, t. In that case we find an explicit solution of the kinetic equation 
using the methods of characteristics. If moreover, the fibre distribution is constant 
in time and space, then we prove the convergence to a parabolic limit. It appears 
to be impossible to prove convergence to the parabolic limit for arbitrary time- and 
space dependent fibre distributions. This confirms numerical observations of Painter 
[21], who investigated the mesenchymal motion model and found interesting cases 
of pattern formation of network type (see Figure 2). In the diffusion limit, however, 
the patterns disappear in the numerical simulation. This indicates that there is a 
significant difference in the asymptotics of the kinetic model and the diffusion limit 
for timely varying tissue networks. 

2. Formulation of the Problem. 

2.1. The Model. We briefly recall the kinetic model for mesenchymal motion from 
[14] for the undirected case. The distribution p{x, t, v) describes the cell density at 
time t > 0, location x G IR" and velocity v G V. Throughout the paper we assume 
that T^ is a product V = [si,S2] x S"~^, where < Si < S2 < oo is the range of 
possible speeds. If si = S2 then we assume si > 0. The fibre network is described 
by the distribution q{x,t,9) with e §"~^, the (n — l)-dimensional unit sphere in 
M" . A schematic of the model is given in Figure 1 . 

fibres - _ \ \ — -^ 



directional change 

- movement direction 
degraded fibres 

Figure 1. Schematic of the model (1) for cell movement in net- 
work tissues, including directional changes, contact guidance and 
fibre degradation. 

The model for mesenchymal motion from [14] reads 
dp{x,t,v) 



dt 



+ V ■ Vp{x, t, v) — —j-ip{x, t, v) + fip{x, t) q{x, t, v), 



^'^ g^^'^^ = K(n„(p(a;, t, v)) - A„(p(x, t, v),q{x, t, 0)))p{x, t)q{x, t, 0), ^^^ 
p{x,Q) ^Po{x), q{x,Q) ^qn{x), 



4 T. HILLEN, P. HINOW AND Z. WANG 

where fj, and k are positive constants. The transport term v ■ Vp indicates that cells 
move with their velocity. The right hand side of the first equation describes the 
reorientation of the cells in the field of fibres. Turning away from their old direction 
at rate /i, they turn into a new direction with a probability that corresponds to the 
fibre distribution q. The new speed is chosen from the interval [si,S2]. The cells 
degrade (at rate k) those fibres that they meet at an approximately right angle while 
they leave fibres that are parallel to their own orientation unchanged. The exact 
definitions of the corresponding terms in system (f) requires some mathematical 
details that are given in the next section. The expressions p, q, n„(p) and Au{p, q) 
are defined in equations (2), (3), (5) and (6), respectively. 

Painter showed in [2f] that the second equation of (f ) arises if instead of ECM 
degradation one assumes that the cells realign the tissue. This would be the case 
for fibroblasts, who do remodel the fibre newtork without destroying it. In that 
case the term —A^ measures the fibre degradation while n„ describes the fibre 
production such that the total amount of fibre mass is preserved. 

2.2. Spaces and Operators. We show in Section 4 that Dirac measures occur as 
meaningful steady states. Hence we need to construct a solution framework that 
allows for measure-valued solutions. Let il = M" be the spatial domain in which 
particles are able to move. 

Let B{V) denote the space of regular signed real-valued (finite) Borel measures 
on V. For p S B{V) let p = p^ — p~ be its Hahn- Jordan decomposition and 
\p\ = p+ + p^ its variation [6]. When equipped with the total variation norm (the 
following notations are used interchangeably throughout the paper) 



\\p\\b(v) = \p\iV) ^ / d\p\iv)^ / H(d«), 
Jv Jv 

B{V) is a Banach space [6, Proposition 4.1.7]. Analogously, S(§"~"'^) will denote 
the Banach space of regular signed Borel measures on S'^^^ equipped with the 
total variation norm. Naturally, we are interested in solutions taking values among 
non-negative measures only. Let 

Xi ==Li(K",B(T/)), 
X2 =L°°(M",S(S""i)), 
X = Xi XX2, 

equipped with norms 

\\p\Wi = / lb(a;)||B(y)da;, 
\\q\\x2 =esssup||g(a;)||e(S"-i), 

x£R" 

ll(p,g)llx-lblk + IMk. 

We denote the positive cones of the spaces Xi, X2 and X by X^^, Xj and X+, 
respectively. We will write 

IIpIIoo =esssup||p(2:)||B(y) 

kSR" 

for those p G Xi for which the essential supremum is finite. 
We define the following operators 
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• The spatial mass density of a velocity distribution, 

-■.B{V)^R, P = p{V). (2) 

Clearly, the operator ~ is Lipschitz continuous. 

• The lifting of a measure on Ei^~^ to a measure on V , 

~:B{W'-^)^B{V), q^m(g)q (3) 

where m is a, probability measure on [si, S2]. If Si — S2, then ^ just maps a 
measure on S"~^ to the same measure on {si} x S"~^. In the paper [14] it was 
taken to be the normalized Lebesgue measure on [si, S2], which corresponds to 
the weight parameter uj defined in [14, equation (4)]. The choice m([si, 82]) = 
1 guarantees that 

l|g||L~(R",e(y)) = Ikllx2- 

In particular, a function that takes values among the probability measures on 
S(§"^^)+ is mapped to a function taking values among probability measures 
on B{V)^. Since " is a linear operator it is Lipschitz continuous. Additionally, 
we use the lifting to connect the measures on V and on S"~^ in a natural way 
as 

dv = m(g) d9, (4) 

• The mean projection operator (for undirected fibres) 



^u{pm - -- , 

PJv 



dp{v). (5) 



For sake of simpler notation and to avoid difficulties when p = 0, we introduce 
the operator 

A:Xi^Li(M",L°°(§"-i)), A(p)=pn„(p). 

Notice that A is linear and if ||p||oo < 00 then 

For sake of completeness we also state the directed version of the operator A, 

A,{p){e)^ f e-/-dp{v). 
Jv m\ 

As said above, existence of solutions is shown completely analogously in the 
two cases. 
• The relative alignment operator again, using the notation from [14] 

A^{p,q)^ [ n„(p)(0)dg(0). (6) 

Similarly to the introduction of A, we will work with 

B:X^Li(M",M), B{p,q) ^ pA^ip,q). 

Notice that B is bilinear and if ||p||oo < 00, then 

ll-B(p,g)||Loo(R„R) < ||p||L-(R",R)||g||x2- 
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The operators A and B are Lipschitz continuous on bounded subsets. 

Let fj, > denote the turning rate and k > denote the rate of fibre degradation. 
The model (f ) can be written as equahty of measures 

dp ^ 

— +vVp^ -tip + ^ipq, 

^^=>,{A{p)~B{p,q))q, (^) 

p{x,0)=po{x), q{x,0) ^qo{x). 

3. Existence Results. To provide a framework for local and global existence of 
solutions we define 

DiA)^{ip,q)eX : VpeXl'}, 
p\_f-v-\7 0\fp\ (8) 



'^' 1 J V ^ J \ 1 

Here V = V^ is interpreted in the sense of weak derivatives of Banach space- 
valued functions. We write / = VxP for a function / € X" if for all test functions 
(j)eW^-^{R",C{V)) 

f{x)-Vx(j){x)dx^ [ pixmx)dxeBiV), 

where the integrals are Bochner integrals taking values in B{V). Observe that the 
domain D{A) is dense in X, as it contains the space C°° {M." , B{V)) x X2 of infinitely 
differentiable functions, which is dense in X [19, Theorem 2.16]. The operator A 
with domain D{A) is the generator of a positive Co-semigroup U{t) on the Banach 
space X (see also Theorem 1 in [2]). 

Notice that the operator —v ■ V is the collisionless transport operator occurring 
in the linear Boltzmann equation which has been studied by many authors, see 
[13, 10], [18, Chapter 13] and the references therein. It generates a semigroup (in 
fact, a group) lAi on the space Xi via 

Ui{t)p(^{x,A)^PQ{x- At,A) -.^ I pQ{x-tdv,dv), (9) 

J A 

for Borel sets A C F. Clearly, the positive cone X^ is invariant under lAi. The 
group lAi preserves the L^-norm while for ||po||oo < c» we have 

\Ui{t)pa\{x,-)^ [ dpo{x-tdv,dv)<{l + ts2\S'''-^\)\\po\\oo. 
Jv 

We denote the semigroup on X generated by the operator A from equation (8) by 
{U{t))t>Q. It has a diagonal structure 

where / denotes the identity on X2. In the operator norm, U satisfies | \U{t) \\c{x) ^ 1 
and for ||mo||oo < 00 we obtain 

||i^(i)uo||co<(l+t.S2|§"-l|)||uo||oo. (11) 

For a pair u = (p, (7) G X define the map 1 1 • | |oo : X — )• [0, 00] by 

Ij-ulloo = IIpIIoo + llgjloo, 
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and set 

I? = {u e X : ||m||oo < oo}. 
For every r > the set 

Dr = {w € X : ||u||oo < '"} 
is closed in X (in particular, the projection of V onto Xi is closed with respect to 
the L^-norm on Xi). Indeed, let Pn € Xi be a Cauchy sequence with ||p„||oo < r 
for all n. Since Xi is complete, it has a hmit p. We claim that ||p||oo < r. Suppose 
that this were not the case, then there would be an e > and a set ^ C M" with 
Lebesgue measure \A\ > such that ||p(a;)||e(y) > r + e for all x € A. But then 
clearly the i^-norm would satisfy ||p„ — p||xi > s\A\ > 0, which is a contradiction. 
Problem (7) can now be written as an abstract Cauchy problem 

u' = Au + F(u), 

u(0) = uo, 

with u = {p, <j), uq = (po, qo) e T). 

Definition 3.1. [20] Let uq — {po,qo) ^ 2?. We say that a function {p,q) ~ u E 
C([0, oo),I?) is a global mild solution if F{u{-)) is continuous and it satisfies the 
integral equation 

u{t)=U{t)uo+ / U{t- s)F{u{s))ds, (13) 

Jo 
where U{t) is the semigroup defined in equation (10). We call a function u = (p, q) : 
[0, T) — >■ I? a classical solution if it satisfies the following properties 

(i) u e Ci((0, T), X) n C([0, T),D{A)), and 
(ii) equation (12) holds. 

Our first result is 

Theorem 3.2. Assume that qo{x,S'^~^) = 1 for almost every x S M", then the 
problem (12) has a unique global positive mild solution for every Uq £ V D X+. 

3.1. Proof of Theorem 3.2. The proof of Theorem 3.2 is established in the fol- 
lowing Lemmas. 

Lemma 3.3. The right hand side of equation (7) defines a nonlinear map 
F : V -^ V, which maps V into itself 

Fiip-.q) \ _ f -t^p + f^pq 



F(P^l) - y F2ip,q) J V 4HP)-B{p,q))q 
The map F is Lipschitz continuous on bounded subsets ofT>. 

Proof. Observe that for (p, q) (z T> the product pq is well defined and 

l|pg||xi < IIpIIli(R'.,r)II<7||x2 = Ibllxilk||x2, 

\\pq\\L°°{R",B{V)) ^ lbllL°°(R",e(y))lkl|x2: 

in particular, 

\\Fi{p,q)\\L--{R'-,B(v)) < 2Ai|bl|L°o(R",e(y))lkllx2- 

For functions (p £ L°°(S"~^) and measures q E S(§"^^) we define the product 
ipq e e(§"-i) by way of 

i^q){M)^ f ^{e)dqie), (14) 

Jm 
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where M C S"^^ is a Borel set. This multiplication extends to functions in L°° (M" x 
S"-i) and L°°(M",6(S"-i)) and we have 

\\m\W2 < II<pIIl-(R"xs— i)lkllx2- 

With ip{e) = K{p){e) - B{p, q) we obtain 

\\F2{p,q)\\oc = ||(A(p) - B{p,q))q\\x, < | |p| U~(R",r)(1 + I kl UJHgllx^, 

showing that F2 takes values in X2. Computations similar to those just carried out 
give the local Lipschitz continuity of F on bounded subsets of V. For example, for 
{Pi,qi), (^2,92) G 2? and IIpiIIxi + Ib2||xi + Ikillxs + Ik2||x2 < K there exists a 
constant C{K) > such that 

llpigi -P2q2\\xi < \\pi{qi ~ g2)||xi + ||(pi -^2)92 1 |xi 
< C{\\pi -P2IIX1 + Iki - 9211x2)- 

We omit the remaining calculations. D 

Lemma 3.4. Equation (12) has a unique local mild solution that remains positive 
for uq G X+ . 

Proof. We set up a Banach's Fixed Point argument, but we cannot work on V 
directly since that set is not complete. Hence we work with Dr for some R large 
enough. For given uq £ V and fixed _R, T > we define 

Er^t ^{ue C{[0,TIVr) : u{0) = uo}. 

This set Er^t is a complete metric space, with the metric given by 

d{u,v) = sup \\u(t) - v{t)\\x- 

te[o,T] 

For a function u e Ert we define 

gu{t) = U{t)uo + f U{t- s)F{u{s)) ds, (15) 

this is again an element of C{[0,T],V) with gu{0) — ug (since V is invariant under 
both the semigroup U and the nonlinearity F). We have for u, v E Er^t 

\\gu{t) - gv{t)\\^ < f \\U{t-s)\\ci^^\\F{u{s))-F{v{s))\\j,ds 
Jo 

< Ct sup \\u{s) ~ v{s)\\x < Ctd{u,v) 
se[o,t] 

(where C is the Lipschitz constant of F) , hence 

d{gu,gv) <CTd{u,v), 

and by choosing T sufficiently small, it can be achieved that ^ is a contraction on 
the space Er^t- If u G "Dr, then we have (see the proof of Lemma 3.3) 

\\F,iv)\\oo<2^iR^ \\F2{v)\U<R\l + R), and 

\\Fiv)\\oo<R^{l + R + 2ti). 

Let u G Er^t- We can estimate equation (15) 

\\gu{t)\\^<{l + ts2\S"-^\)\\uo\\oo+til + ts2\S''-'\)R^{l + R + 2^i). 
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By choosing R > 2||uo||oo (large) and T small, namely 

1 i?-2||uo||c 



T < min 



S2|§"-i|' 2R2{l + R + 2n)j ' 
we can achieve that 

sup \\gu{t)\\oo < \\gu{T)\\^ <R, 

te[o,T] 

for all u G Eji^t- Hence the contraction Q maps the complete metric space Eh^t 
into itself and hence has a unique fixed point by the Banach Fixed Point Theorem. 
The positivity of solutions follows from the fact that the nonlinearity F is of 
multiplicative type. If either p or q becomes zero on a set at some time, the left 
hand side of equation (12) is non-negative. D 

Concerning the global existence of solutions, if Tmax{uo) < oo, then by [22] 

lim ||u(i)||x = oo. 

However, in our system (7) a blow-up in finite time cannot occur as the following 
lemma shows. 

Lemma 3.5. Let {p,q){t) be a mild solution of equation (7) (equivalently, of (12)) 
taking values in V D X+. Then for all t € [0,Tmax) o,i^d almost every x S M" we 
have 

g(x,i,§"-i) = l, 
and there exists a constant C > such that 

lb(i)lk = |bo||xi, and |b(i)||oo<|"^^||bo||ooe^*. 

Proof. Let (p, q) be a mild solution of equation (7) taking values in X+. The 
second component of equation (13) reads 

q{x, t) = Iqo{x) + / Ik{A{p) - B{p, q)q) ds, 
Jo 

where / denotes the identity. We evaluate this relation at §"~^ and use the fact 
that 

k{A{p) - B{p, g))9(§"-i) = / A(p)(0) Aq{x, 9) - B{p, g)g(§"-i) 

= i?(p,g)(l-g(§"-i)). 



We obtain 

1 - g(a;,t,S""') == 1 - qoix,S"-') - k I B{p,q){l - q{x,s,S"~')ds 

Jo 

We apply Gronwall's lemma and obtain 

/■* 



l-g(x,i,S"-0 = (l-go(a;,S""0)exp(-^y i?(p(x, s), g(x, s)) ds 1 . 

The integrand is positive and bounded, hence by the assumption on qo we get 

g(a;,t,§"-i)==l (16) 

for almost all x E M". For the L^-norm of p wc notice first that since p is positive, 
it satisfies 

p{x,t,V) = \\p{x,t)\\B{V)- 



10 T. HILLEN, P. HINOW AND Z. WANG 



We evaluate the first equation of (7) and obtain as a consequence of (16) 
d_ 

Integrating tliis equation over M" gives 
d 



^p(x,t)+V- ij vAp{x,t,v)\ =0. 



, , p{x, t) dx = — V • / u dp(x, t,v) \ dx = 

by the divergence theorem. For the L°°-part of p we use that Fi(0, g) = and the 
following fact 

\\Fl{pi,q) -^l(P2,g)||L-(R",B(y)) < Ai(l + I|9||X2)I|P1 -P2||L°o(R",f3(y))- 

We estimate from (13) 

lb(i)l|oo<|l/|(^lbo||oo+^ ||Fi(p(s),g(s))|Uds 

= l^l(lbo||oo+y I |Fi(p(s),g(s))-Fi(0,9(s)) Hoods') 

<|V^l(lbol|oo+A*(l+ sup ||g(s)||xj / Ib(s)lloods). 
\ se[o,t] Jo y 

This inequality warrants application of Gronwall's lemma 

Mt)\\^<\V\\\p,\Ue^' 

with a suitably chosen constant C . By the density of the domain £'(.4) in I? n X 
and because of the continuous dependence of the solution on the initial datum we 
obtain the desired estimates for arbitrary initial data [pq, go) S T^ C\ X+. D 

Combining Lemmas 3.4-3.5 we conclude the proof of Theorem 3.2. 

4. Steady States. In numerical simulations by Painter [21], shown in Figure 2, we 
find interesting network patterns which form from random initial data. Numerically, 
these patterns do not change after they have been established. We expect that the 
system (7) allows for these network patterns as steady states. In this section we will 
develop a theory of pointwise steady states which are candidates for the observed 
network patterns. 

To describe steady states of (7) we introduce the bilinear turning operator 

C ■.B{§"-')xB{V)^B{V), C[q]{p)^qp-p. 

Observe that in contrast to the paper of Chalub et al. [3] the turning kernel does 
not depend explicitly on w', i.e., the cells are reoriented regardless of their original 
orientation. For p G kei C[q], we have 

p^qp. 

Hence the orientation of the cells in a steady state is entirely given by the fibre 
distribution q. This reflects the fact that a perfect alignment of the cells with the 
underlying fibre network and only such a perfect alignment remains invariant under 
the turning operator C. 

The trivial steady state is a uniform distribution of fibres and cells: 
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Figure 2. Typical simulation of network formation for model (7). 
The left figure shows the overall cell density p(x, €) at a time where 
the steady state has almost been reached. Light (red) color indi- 
cates high cell density and dark (blue) color indicates low cell den- 
sity. The figure on the right shows the underlying network, where 
the small bars indicate the mean direction and the gray color de- 
scribes the degree of alignment. Light gray indicates highly aligned 
tissue, whereas dark gray/black indicates close to uniform distri- 
bution of directions. The simulations were done by K. Painter, and 
are described in detail in [21]. We are grateful to K. Painter who 
allowed us to use this figure for illustrative purposes. 



Lemma 4.1. (Homogeneous tissue) For every constant g > the pair 



q{x) 



AO 



|§»- 



p{x) ^ gq^ g 



|S"-i| 



dv 



0^ 



is a steady state of (7) in L°°(M",B(y) x B{E"-'^)). The only steady state of this 
type m I? n X is obtained for g = 0. 

Proof. If q = d0/|§"~^ I and p — gq, then p = g and p = pq. The right hand side 
of the first equation of (7) is zero. For the second equation, we need to compute 
A(p) — B{p, q). We have 



A{p){x,0)^g 



V 



V 

M 



dm = /3 



for a /3 > 0, which is independent of 9 and x. We obtain 

d9 



B{p,q){x) 



/3t 



/?• 



S"-i 



Hence A{p) — B{p,q) — and the right hand side of the second equation of (7) 
is zero as well. Notice that the measures Av and A9 are coupled in a natural way 
through (4). D 

To find other steady states, we need a weak formulation. 

Definition 4.2. We say that (p, g) G PnX is a weak steady state of (7), if for each 
pair of test functions 

cf) e M^i'i(M", C{V)), ^ e Co(M", C(S"-i)), 
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(where Cq denotes functions vanishing at cx)) we have 

V ■ 'S/4>{x, v)p{x, dv) dx = 
V 

X, v) {—^p{x, Av) + fip{x) q{x, dv)) da;, (17) 

V 

iA{p){0) - Bip, q)) t^ix, 9)q{x, d0) dx = 0. (18) 

-1 

Notice that in this definition, and tp are real-valued functions in the variables 
V and 0, respectively, hence the integrals on V and §"~^ make sense. In the next 
Lemma we study the biologically meaningful case of a network completely aligned 
in a single direction. 

Lemma 4.3. (Strictly aligned tissue) Assume a preferred direction 7 S S"~^ is 
given and g > is a constant. Let S^ denote the Dirac mass on S"~^ concentrated 
at 7. Then 

p{x) = gq, q{x) == '^ ^ ~^ 

is a weak steady state m i°°(M",6(y) x B{E>"-^)). 

Proof. Since p e ker>C and since it is spatially homogeneous, equation (17) is 
satisfied. To study (18) we first compute the following integrals on §"^^ 



'V 



iS^+d_^)idv)^g\0-j\, 



A{p){x, 0)^{x, 0) dq{0) = |(V(a:, 7) + ^{x, -7)) 

B{p,q){x) = ^ f \0-j\{6,+6.^){d0)^g, 



and 

B{p,q){x) f „,,i^a\ix ,x ^/J^^ _ ^ 

2 
Hence 

(A(p)(x, 0) - B[p, q){x)) V;(x, 0) dq{e) - 



V'(x, 0) {6^ + 5.^){d0) = ^(V'(x, 7) + i^{x, -7))- 



for all xeW\ D 

4.1. Pointwise Steady States. In the preceding Lemmas we identified two simple 
homogeneous steady states. A full analysis of other steady states at this level is 
difficult, since the very weak formulation of measure valued solutions allows too 
many degrees of freedom. We rather specialize to the study of pointwise steady 
states as defined below. With pointwise steady states, we can combine the previous 
two Lemmas and design networks of aligned tissue with patches of uniform tissue. 
A schematic of the steady states which we construct here is shown in Figure 3. 

Definition 4.4. We say that (p, g) G 2? n X is a pointwise steady state of (7), if 

1. [p, q) is a weak steady state. 

2. p{x), q{x) is well defined for each x € E". 

3. For each test function ^ e C(S""i) and each x e M" 

{A{p){0) - B{p, q)) ^{0)q{x, d0) = 0. (19) 



CELL MOVEMENT IN NETWORK TISSUES 



13 




Figure 3. Schematic of steady states. Figure A and B show typi- 
cal fibre alignment and cell alignment for the homogeneous steady 
state (in A) and the strictly aligned steady state (in B). Figures 
C and D show the geometric construction that underlies pointwise 
steady states. The fibres and cells are aligned tangentially along 
the curves ct^ and uniformly inside the domains fli . Figure C shows 
a pattern without intersections, while Figure D shows intersections. 
One of the intersections has been blown up to illustrate how a three 
pointed star of 120° angles can arise (see Corollary 1). 



4. For each test function <I> € C{V) and each x E 



$(?;) {—^p{x, dv) + ^.p{x) q{x, dv)) — 0. 



(20) 



Remark 1. : (a) An immediate consequence of item 1. and 4. in this definition 
is that pointwise steady states satisfy 



V ■ \'(j){x, v)p{x, dv) dx = 0, 



(21) 



for each test function (j) £ W^^'^iW\C{V)). 
: (b) Another immediate observation is that the homogeneous steady state from 
Lemma 4.1 and the completely aligned steady state from Lemma 4.3 are point- 
wise steady states. 

In the following we classify pointwise steady states in M.'^. It turns out that the 
above definition allows for patchy steady states and for steady states of network 
type. Patchy steady states include patches of uniform tissue surrounded by areas 
of strictly aligned tissue, i.e. a combination of the above two types. Network type 
steady states arise if the areas of aligned tissue form a connected network of curves 
with intersections and branches. For network type steady states, we will classify 
possible intersections of network fibres. We are able to explicitly treat intersections 
of up to four directions and we find a general algebraic condition for networks with 
intersections of higher order. 
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4.2. Patchy Steady States. Assume a set of smooth curves (7^, i = 1,...,N 
separate M^ into disjomt open sets fli,i = 1, . . . , fc. Assume these curves ai have 
finite length, no intersections but they might be closed. Assume p and q are uniform 
inside each patch 

P{^)^P^Tyi, ^(^) " |§i| ' lixen^, i = l,...,fc (22) 

with Pi = if |0i| = oo. Since we are interested in pointwise steady states, we 
need to define {p{x), q{x)) for x E [J- oi. For each i = 1, . . . , iV we denote the unit 
tangent vector at a; G ct^ by 7i(a;), where we will suppress the argument x whenever 
possible. We define 

9i(a;) = -^{^-^dx) + ^liix))-, P{x) = gtq{x), for x e a,. (23) 

and Qi > 0. 

Lemma 4.5. The weak steady state defined by (22) and (23) is a pointwise steady 
state. 

Proof. {p{x),q{x)) are defined for all a; G M^ and as shown in the proofs of 
Lemma 4.1 and 4.3 the conditions (19) and (20) are satisfied for all x E M^. We 
only need to show that (p, q) as defined above is a weak steady state, i.e., we need 
to confirm condition (21). We find 

V ■ V(/)(a;, v)p{x, dv) dx 

2 JV 

k k 

— /^ / 4'{x,v)v ■ S/ p{x, dv)dx — y^ / n ■ v(t>{x,v)p{x, dv)da, 
j^i Jn, JV j^;^ Ju, Jv 

= 

The first integral vanishes since VxP{x, dv) — QinVti. The boundary integrals are 
zero, since we assumed that on cr^ the fibre orientation is tangential, i.e. n ■ v = Q 
for all V G supp p{x, i), where n denotes the outer normal of cr^ at a; G cr^ . D 

The above lemma allows for patches of uniform tissue surrounded by aligned 
tissue. These could also be called encapsulations, as seen for many tumours in 
tissue. A schematic of patchy steady states is given in Figure 3C. The steady states 
in Lemma 4.5 do, however, not allow for intersections of the curves ai so that they 
become of network type. To obtain network steady states, we need to study possible 
intersections in more detail. 

4.3. Symmetric Intersections. To study multiple directions we introduce two 
abbreviations. For a given vector 7 G S^ and a real valued function ^ on S^ we 
define the notation 

<5h :- \{5-, + 5,), ^{yi\) := ^(^(-7) + ^^(7)). 

We first consider the intersection of two directions 71,72 G S""^ with 71 ^ ±72 
and with different weight a G (0, 1). For a given a; G M^ we define 

q(a:) :=a(5|^^l + (1 -a)(5|^,| and p{x) ^ Qq{x), (24) 

where we set g = 1 without restriction. 
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Lemma 4.6. Assume {p, q) is a weak steady state of (7) and at x it is of the form 
(24)- It can only be a pointwise steady state, if the directions 71 and 72 have equal 
weight, i.e. if a — ^. 

Proof. We only need to check condition (19) of Definition 4.4 at the intersection 
point X. For this choice of p and q we find 



Hp) 



V 

M 



dq ^ a\d ■ ji\ + {1 - a)\0 ■ -f2\ (25) 



and 



B{p,q)^ / A{pmdq 

= a^|7i7i| + 2a(l - a)|7i72| + (1 - a)^|7272| 

= 2a2 -2q:+ l + 2a(l-Q:)|7i72| (26) 

where we used the fact that |7i7i| = l,i = 1,2. To check condition (19), we need 
to test with a test function ^ e C(§^): 

A{p){e)^{0) dq = a{a + (1 - a)|7i72|)*(l7i|) 

SI 

+ (1 - a)(a|7i72| + 1 - a)^(|72|), 
B{p,q)^{9)dq^B{p,q){a^{\j,\) + {l-a)^{h2\))- 

Hence to satisfy (19) for any test function, we need to satisfy 

a + (1 - a)|7i72| = B{p,q) = a|7i72| + 1 - a. (27) 

Comparing the first and last term, we obtain the equation 

2a-l = (2a-l)|7i72| 

which is satisfied only if a = i. Notice that we assume I7172I 7^ 1- For a = ^ we 
find 

B{p,q) = - + ^\lil2\ 

and hence the condition (27) is satisfied. D 

Next we study the general case where at a given point a; G M'^ we have an 
intersection of A^-different directions 71, . . . ,7jv G S^- We study N directions with 
equal weight: 

9(2;) := ^(^|7i| +--- + '5|^„|) and p{x)=q{x). (28) 

To decide if this intersection can be a pointwise steady state, we define a matrix of 
pairwise projections: 

T■■^{\l^l^)^,l = l,...,N■ (29) 

Theorem 4.7. Assume {p,q) is a weak steady state of (7) and at x it is of the 
form (28). It can only be a pointwise steady state, if the corresponding projection 
matrix T has an eigenvector (1, . . . , 1)"^. 

Proof. Again, we only need to check condition (19) of definition (4.4) at the 
intersection point. For the above choice of p and q we find 

1 ^ 

Hp) ^ ^T.\(^l^\ (30) 
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and 

1 ^ 
B{P-^^)-^2Y.\^^^^\- (31) 

ij = l 

Applied to a test function 4* e C(§^) we obtain 

£ A(p)(0)*(0)dq= ^f] fEl7,7d*(l7,l)j (32) 

and 

/ B{p,q)-^{e)dq = -^ Y. l7nd ^E*(l^^l)- (33) 

•^^^ ij=l fc=l 

To satisfy condition (19) the right hand sides of (32) and (33) have to coincide for 
any test function. In particular we need to satisfy 

N N 

Y. \^i^-\ ^N^ l^^"^^' I' ^°' "'"''^ l = l,...,N (34) 

z— 1 ijj — 1 

This condition implies 

N N 

X! \^i-^i\ = X! I'^'^'^*! 1°"^ '^^^^ l,k = l,...,N. (35) 

1=1 4=1 

It can be directly verified that condition (35) implies (34) and it also implies 
(32)=(33). Hence (35) is the limiting condition. This condition implies that the 
row-sums of the matrix T are all identical, and since F is a symmetric matrix, the 
column sums also have the same value. In other words (35) is equivalent with the 
statement that F has an eigenvector (1, . . . , 1)^. D 

A schematic of a steady state with intersections is shown in Figure 3D. 



Remark 2. Notice that a related matrix to F is well known in linear algebra: the 
Gram matrix 

plays a role in coordinate transformations and the square root of the Gram deter- 
minant is a measure for the volume element spanned by the vectors 71, ... , 77V. 

Example 4.8. As an example, we apply this general result to the two-directional 
case studied in Lemma 4.6. For two directions we have 

^^(hn,! 'T'). - r(;)^a+i,..i,(; 

For three directions we obtain an interesting result: 

Corollary 1. Assume (p, q) is a weak steady state of (7) and at x it is of the form 
(28) with N — 3. It can only be a pointwise steady state, if the three directions have 
equal angle, i.e. I7172I = I7273I = l737i|- 

Proof. We use the criterion from Theorem 4.7. The vector (1, 1, 1)"^ is an eigen- 
vector of F if 

1 + I7172I + I7173I = 1 + I7172I + I7273I = 1 + I7173I + I7273I 
which implies 

I7172I = I7273I = l737i|- 
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D 
A three pointed intersection has been illustrated in Figure 3D. 

The classification of intersections of four directions is a bit more complex. 

Corollary 2. Assume (p, q) is a weak steady state of (7) and at x it is of the form 
(28) with N — 4. It can only be a pointwise steady state, if the pairwise equal angle 
condition (36) is satisfied. 

Proof. To illustrate this case we introduce another abbreviation 

Qtj ■■= |7j7jI- 
The corresponding projection matrix for four directions reads 

/ 1 512 ffl3 514 ^ 

p _ 5l2 1 .923 524 

5l3 523 1 534 

V 514 524 534 1 J 

and the eigenvalue condition is given by 

1 + 5l2 + 5l3 + 5l4 = 1 + 5l2 + .923 + .924 = 1 + 5l3 + 523 + .934 
= 1+514 +524 +534- 

Hence we obtain six unknowns and three equations, which will not give us such 
a complete solution as for three directions. We can, however, reduce the above 
condition to a set of pairwise equal angle conditions 

I7172I = I7374I, I7173I = I7274I, I7174I = |7273|- (36) 

If all angles are 7r/2 then this condition is satisfied. D 

4.4. Unsymmetrical Intersections. In the numerical simulations shown in Fig- 
ure 2 we observe intersections that are asymmetric in the sense that for a direction 
7i the opposite direction —71 is not seen. Indeed, at an intersection point x various 
fibres come together. This means that at this point x the network has a number 
of directions 71, ... , 7Ar, but the opposite directions are missing (it could of course 
happen by chance that 7^ — —ji for some i, j). Even though we assume that the dis- 
tribution function q is symmetric almost everywhere, we find exceptional points at 
those intersection points. In this section we show that unsymmetrical intersections 
can arise as steady states in the framework developed here. 

Assume at a given point a; G K^ we have an unsymmetrical intersection of A'^- 
diffcrent directions 71, ... , 7Ar G §^ with equal weight: 

1 

To decide if unsymmetrical intersections can arise as pointwise steady states, we 
carry out the same computations as in the previous section. It turns out that the 
computations change only marginally and we omit the details here. For example 
formulas (32) and (33) use ^(74) instead of ^(|7i|). This implies that the conditions 
for their existence remain the same. We summarize: 

Theorem 4.9. 1. Assume {p, q) is a weak steady state of (7) and at x it is of 
the form (37). It can only be a pointwise steady state, if the corresponding 
projection matrix T has an eigenvector (1, . . . , 1)^. 
2. If N — 3 then (p,q) can only be a pointwise steady state, if the three directions 
have equal angle, i.e. I7172I = I7273I = l737i|- 



?(^) ■= l^('^7i ^ ^^7iv) and p{x)^q{x). (37) 
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3. IJ N = A it can only be a pointwise steady state, if the pairwise equal angle 
condition (36) is satisfied. 

Remark 3. : (a) Indeed, unsymmetrical intersections of three directions with 
angles of 120° seem to be typical building blocks for the network shown in 
Figure 2. 
: (b) Although other intersections (symmetric or asymmetric) do exist theoreti- 
cally, they are rarely seen in simulations. This raises the question of stability 
of these steady states. We defer this question to future studies. 

4.5. Other Steady States. We consider more general steady states where the cell 
distribution is a multiple of the lifted fibre distribution, that is 

p{x) = g{x)q{x), 

where g g i°°(M") is the density of cells (or even q & L^ f^ L°°(M")). The minimal 
condition for such a pair to be a steady state is 

K{p){e)=B{p,q), 

in particular, the left hand side is actually independent of 0. Because of the linearity 
of the operators A and B, this condition becomes 

Q{x)K(q{x)m = Q{x)B{q(x),q{x)). 

Wherever g 7^ this condition can be stated as 



V 



g(a;,dw) 



q{x,Av)q{x,de), (38) 



for almost all x e M". Now let us try the following ansatz 

q{x) = /(x)(5|^(^)| + (1 - /(a;))!], 

where < f{x) < 1 and E is the normalized Haar measure on S"~^. Notice that 
even the predominant direction 7 may depend on x at this point. A calculation 
gives 

V 



q{x,dv) = / \e-'ip\q{x,dip) 



l^ll 



The last term in the second line is independent of 6 because of the rotational 
invariance of S. Hence we have 



fix) + (1 - fix)) / \0-^\ S(dV') =: C. 



Cqix,de)^C qix,Ae) = C. 

Observe that 

Qix)qix) = gix) 

and the right hand side of the first equation of (7) vanishes. Hence the condition 
for Q is 

Qix)Vq{x) + qix)VQix) = 0, 
for almost all x G E", which can be written as 

S/iQix)qix)) = 0, (39) 

In summary, if we have found q that satisfies (38), then g can be determined from 
the differential equation (39). 
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5. Analysis For Given fibre Distribution. In this section we assume that cells 
do not remodel the fibre network and that q{x,t) is a given distribution. The p- 
equation of system (7) has a simple structure for given q. For classical solutions we 
can use the method of characteristics to find an explicit solution. For given v € V, 
the characteristic equation is ^x{t) = v. Hence the characteristic through xq S M" 
is given by x{t) = xq + vt. We can write the first equation of (7) as follows 

j^p{x{t),t) + M^it), t) = mi^it), t)p{x{t),t), (40) 

We evaluate equation (40) at V and obtain 
-p(a;(i), t) = j^P{xit),t, V) = ^Mx{t),t, V) + fip{x{t),t)q{x{t),t, V) = 0, 

where we have used the fact that q{x{t), t, V) = 1. Hence p{x{t),t) is constant along 
cliaracteristics. Equation (40) is equivalent to the equation 

e~''*-{p{x{t),t)e''*) = fiq{x{t),t)p{x{t),t). 
Integrating the above equation with respect to time, we obtain 

p{x{t),t) ^ e-^V(a;o, 0) + fie-^^'pixit), t) f e''"g(x(s), s) ds. (41) 

Jo 

For a given {x,t) e M" x IR+, we find the anchor-point xo{v) ^ x — vt and the 

corresponding backward characteristic in the direction v 

x{s) — X — vt + vs. 

Applying this to equation (41), we have 

p{x,t) = €''''' Poix-tdv) + fip{x,t) e-''''*-'^q{x- {t- s)dv,s)ds. (42) 

Jo 

This is an equality in the Banach space B{V) and the term po{x — tdv) on the 

right hand side has to be interpreted as the v-shifted measure defined in equation 

(9). The same notation apphes to q. We evaluate this measure p{x,t) at V, i.e., we 

compute p{x, t) = p(x, t, V), and obtain 

p{x, t) = e-^'^paix - Vt, V) + fip{x, t) j e-'''^^-'^{x - V{t -s),s, V) ds, 

this is an equality between real numbers. The measure q is non-negative and for 
fixed w e y we have q{x — w(t — s), s){V) = 1, and 

K{x, t)^fi I e"^(*"")g(a; - V{t - s), s, V) ds > 0. (43) 

Thus we get 

(1 - K{x, t))p{x, t) = e-f^'poix - Vt, V). 
If K{x, t) y^ 1, then we can solve for p as 

P{x,t)^ ^_'''j^^^^f oix-Vt,V). (44) 

Then p can be used in (42) to find an explicit solution 
p{x,t) =e"''*po (x-tdv) 

up->^* /■* , ^ (45) 

+ 1 r^( ,M ^ - ^*' ^) / e-''(*-^)g(x -{t-s) dv, s) ds. 
1-K{x,t) Jo 
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Notice that this solution only depends on the initial condition pq and on the fibre 
distribution q. To clarify, equation (45) is again an equality in B{V) and the right 
hand side is of the type "measure + number • measure" . 

Equation (43) simplifies drastically in the special case of constant fibre distribu- 
tion q{x,t) = q. In this case, it follows that 

ft 







~Kt-s)~ 



q(y) ds = ^J■ I e 




-t^(t~^}ds = l 



-fit 



K{x,t) 

Equation (44) becomes 

p(x,t)^Po{x-Vt,V), (46) 

Equation (46) deserves some interpretation, p is the mass density of particles of 
all velocities at point {x, t), whereas pQ{x~Vt^ V) integrates the initial condition over 
the domain of dependence of the point (x, t), the set {x — tv : v G V}. The velocity 
distribution at (x, t) arises by following all characteristics through (x, t) backwards 
(see Figure 4). We call (46) a generalized Huygens principle. The solution p{x,t) 
from (45) can then be written entirely in terms of the initial condition 

p{x, t) = e-f^'poix -tdv) + {l- e-^*)po(a; - Vt, V)q. (47) 

Using equation (46), the explicit solution can also be written as 

p{x,t) ^ e'-^'^poix - tdv) + {1 - e-''*)p{x,t)q. (48) 

Hence the solution is a convex combination of the initial condition pq and the current 
amount of cells p redistributed with respect to the "controlling" distribution q. We 
will use this observation in the next section to rigorously prove convergence of the 
parabolic limit. 




Figure 4. The domain of dependence of the point x at different 
time points is shown as a thick solid line on the x-axis. In this 
example V = [si, S2] x S"^^ with si > is an annulus. 



It is interesting to understand the biological meaning of these explicit solutions. 
Equation (48) tells us that eventually cells will completely align to the given network 
structure. This has similarities with glioma cell invasion along white matter tracks 
of the brain [17]. Equation (45) has a similar interpretation. In contrast to (48), 
here the fibre distribution varies in time and space. The integral terms in (45) and 
(43) denote a temporal average over the history of the fibre orientation, where the 
influence of the history is exponentially damped. Then (45) can be understood 
as cells that try to align with the tissue while the tissue is changing. We see 
a biological analogy with wound healing, where fibroblasts constantly modify the 
collagen network while immune cells move through this fibre scaffold and heal the 
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wound [7] . Notice that it is not our intention here to model brain tumours or wound 
heahng. These examples are only used as analogies and thought experiments. It 
might be useful to make our model available to these processes in the future. 

5.1. The Parabolic Limit Problem. As shown in [14], we can formally derive a 
diffusion limit equation from equation (7) under suitable scaling of space and time. 
Let X and v denote reference length, and speed, respectively, with the dimensionless 
quantity 

V 

e — — 
being small. We introduce rescaled variables as follows 









r = eh 


* ex 

V 


and V* = 


V 

v' 




This 

and ' 
tion 


gives 
we obtain. 


upon 


d 
dt 

dropping 


dt*' 
the asterisks. 


V 

the reduced parabolically scaled 


equa- 








s^ ;;+... vp. = - 


-fiC[q]{ps), 




(p.: 



Pe{x,0) ^poix) eX'nX+. 

Simultaneously, we consider the limit problem 

e(x, 0) = poix, V) = poix) e L''+{W\R), 
with ||g( • , 0)||oo < oo and where the diffusion tensor is given by 

D[q]^- / v<Sivdq{v). (49) 

^J■Jv 

The formal derivation of the limit problem and the diffusion tensor in equation (49) 
was carried out in [14, section 4], see in particular equations (29) and (41) in that 
paper. We therefore omit these calculations here. Notice that D[q] can be written 
as the scaled variance-covariance matrix Y{q) of q, 

D[q] = - f ' f {s9) «) {sO) dq{9) dm{s) = crV(q), 
M Jsi Js— 1 

where 

^ ' „2 



<J=- s^dm{s), V(g) = / 9 ® 6 dq{e) 

and we have used equation (3). 

We define the notion of a weak solution of equation (Po). 

Definition 5.1. Let T > be given. We say that g e W^-^{[0,T], W'^-^{R'^,R)) is 
a weak solution of (Po) if the following holds 

Q(x,0)(j){x,0)dx — / / g(x,t)^—(x,t)dtdx 
Jr" Jo dt 



/ / D[q]{x,t)yg{x,t)-V(j){x,t)dtdx 

Jr" Jo 
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for all test functions </> € C^{[0, T] x E") with 0( • , T) = 0, and in addition 

g{x,0) =po{x) 

for almost every x G M". 

The tensor D[q\ in equation (49) is positive definite as long as the support suppq 
is not contained in a lower dimensional great sphere. To see this we take a £ M" 
and study 



a^Dlqja = - i {v ■ af dq{v) > 0, 
A* Jv 



provided that suppq is not contained in (a)^ n §" ^ for any a G M". In this case, 
we have the existence of weak solutions [12, 23]. 

5.2. Convergence Result. The parabolic diffusion limit for chemotaxis was rig- 
orously studied by Chalub et al. [3]. It was assumed that there exists a bounded 
equilibrium velocity distribution F{v) € L°°{V) that is independent of space, time 
and the distribution of the chemical signal, [3, Assumption (AO)]. The assumption 
(AO) in [3] corresponds to our assumption (AO) below for the case that the equilib- 
rium distribution of the turning operator is a given function/measure on S"~^ (and 
independent of x and t). The difference arises from the fact that F is uniformly 
bounded while g is a Borel measure. 

Since we are now equipped with a suitable functional analytical setting, we will 
rigorously study the convergence to the parabolic limit. However, as shown numeri- 
cally in Painter (see [21, Figure 9]), the phenomenon of network forming patterns is 
lost in the diffusion limit, hence we do not expect that convergence to the diffusion 
limit is true in general. We assume that q is constant in space and time 

q{x,t) = qeB{§''-^) (AO) 

for all X E M" and t > and that q is symmetric with respect to 9 i-^ —6. 

Theorem 5.2. Let assumption (AO) hold and fix T > 0. Let (p^)^>o be the fam- 
ily of solutions to problem (P^) and g the weak solution to problem (Po) (in the 
sense of Definition 5.1). Then, after possibly extracting a subsequence we have the 
convergence 

Pe -^ gq 

in the weak* topology on the .space L°°([0, T], Xi). 

Proof. Let p^ denote the B{V)-va,\ued solution of equation (P^). We solve this 
equation as we did in Section 5, observing the new scaling with respect to e. After 
dividing equation (P^) by e^ and applying (47), we find 

p,{x,t)^e-^%(x-^-^j +{l-e-7^')po(x~^,Vjq. (50) 

The family {p^)^>Q is uniformly bounded in i°°([0,T],Xi) since ||p(-,i)||xi = 
IIpoIIxi- Hence there exists a weak* -convergent subsequence, say Pe ^ p* as £ — >■ 0. 
Taking the B(V) norm in equation (50) and then taking the supremum over all 
a; e R" and i e [0, T] gives 



\\Pe{x,t)\\B(v) 



< 



tdv 
Po ■ 



B{V) 



Po[x ,V 



<2||po||oo. (51) 
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and hence ||p*(- ,i)||oo < oo. Using equation (46) in the rescaled coordinates we 
rewrite equation (50) as 

Pe{x,t) =e~^* ipoix j -~p^{x,t)q] +p^{x,t)q. 

Sending e to in this equation we see that p^,{x,t) — g{x,t)q for an appropriate 
function g G L^{W^ x [0,T],R) with ||p(- ,i)||oo < oo. It remains to prove that g so 
defined satisfies the paraboUc hmit problem (Pq)- To this end we define a residuum 
r^ and obtain with (50) 

Pe-peQ e-^' / / tdv\ ( Vt \ . 

Observe that f^ = and for e > 

e 
By a similar argument as for p^ in (5f ), we get 

||rs(x,t)||e(y)<2|bo||oo. 

Hence there exists a weak*-convergent subsequence r^ ^ r*. Finally, let (/? S 
Cl{M."- X [0,T],M) be a test function and observe that 

^ / / -^{x,t)ip{x,t) Ax At 
Jo ■/»" ot 

dip 



Pe{x,t)ip{x,t)dx 



Pe{x,t) — {x,t)dxdt. 



Since the right hand side converges to zero as £ — >■ 0, so does the left hand side and 
we have that 

in the distributional sense. 

We divide equation (P^) by e and obtain 

dps ^ 

e— +v-Vpe = -fJ-re- 

Now we let e — ^ 0, divide by /i and we obtain the following representation of the 
limit of the residuum 

r, = vV{gq). (52) 

We evaluate equation (P^) at V and obtain the conservation law 



.2^P£ 

dt 



J vd{er,+p,q)]^0. (53) 



By the symmetry of q, we have 

vdq{v) — 0. 



'V 

We divide equation (53) by e^ and let e — > and obtain 

dg 



dt- -v-y^«dr.(.), 
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where 

dpe _^ dg 

dt ~^ dt 
in the distributional sense. Using the representation (52) we obtain 

^=V-[-/u(8)W dq{v) Vg 

Hence g satisfies the hmit equation (Po)- Q 

6. Discussion. In this paper we consider mathematical properties of a model (1), 
or equivalently (7), that describes mesenchymal cell movement in tissues. The 
model was developed in [14] and has been analyzed from various angles in recent 
papers, [21, 4, 5, 25]. Through the previous analysis it became evident that a 
solution framework is needed which allows for measure valued solutions. Here we 
develop such a framework and prove global existence of solutions in X. We have 
used semigroup methods, since they provide a dynamical systems point of view, and 
we can use this framework for linear stability analysis in future work. Alternative 
methods to show existence include energy methods as developed by DiPerna and 
Lions [9]. 

We were able to find non-trivial measure- valued steady states, which correspond 
to homogeneous distributions, or to aligned tissue, or to patches of uniform tissue 
with a network separating these patches. We found that pointwise steady state 
show network properties as observed numerically. We also found that, although the 
system has been formulated symmetrically, we can have unsymmetrical intersection 
points. This confirms the interpretation in [14], where it was suggested that a 
network made from undirected fibres can have characteristics of a directed network. 
The complete identification of steady states of (7) is an interesting open question. 
Furthermore it would be interesting to see whether solutions of (7) converge to 
steady states or to traveling wave solutions as i — > oo. The existence result in X 
opens the door to a rigorous linear stability analysis of steady states. This endeavor 
is left to future work. 

The convergence to the parabolic limit is a standard feature of kinetic models 
and it has been studied in many publications (see references in the text). Our 
approach here extends known results to measure-valued solutions. Furthermore, we 
formulate an explicit solution which shows that the solution basically is a convex 
combination of initial data and its velocity-mean- value. The mean value then is close 
to the parabolic limit. We also give an argument that the rigorous convergence to 
a diffusion limit might only work for constant tissue. 

Here we did not discuss the biological modelling of (7). We would like, however, 
to discuss the biological assumptions and propose various extensions, which could 
lead to more realistic models. 

One possibility is to introduce birth and death processes for the cells into the 
model. For example, it is known that growth factors can be bound to the fibres 
which promote proliferation. Also, harmful substances can be found in the fibre 
network, possibly killing the cell. To model these effects, a term G{p) would be 
added to the right-hand side of the first equation of (7). If G{p) satisfies certain 
growth bounds, the global existence of solutions will continue to hold. 

A second possibility is to allow diffusion of p with respect to both x and v. 
Cells are likely to undergo some random walk and may also change their velocity 
randomly (a perfect alignment of cell velocities will disperse). Diffusion with respect 



CELL MOVEMENT IN NETWORK TISSUES 25 

to the X variable is easily modeled by adding a term of the type —Dr^A^p to the 
left-hand side of the first equation. Also, diffusion in the velocity can be modelled 
through an additional diffusion term of the form —DyAyp (see also Dickinison [8] 
for chemotaxis) . For these cases we expect a smoothing property of the linear 
semigroup and the totally aligned steady states will no longer exist. 

Another possibility to expand and make the model more realistic would be to 
give the fibres some elasticity and to let the fibres be moved by the cells. Also, the 
cells should chose their new speed not randomly, but according to some "stiffness" 
of the neighborhood they are currently in. For example, a cell that has to cut a 
lot of fibres in its way should slow down, while a cell that is aligned well with the 
network can gain speed. Obviously, these are intuitive ideas, and would have to be 
supported by biological evidence. 

In model (7) we implicitly assume that the protease is released locally at the lead- 
ing edge of the cell. In the literature, however, various protease cutting mechanisms 
are discussed [11] and more detail of the cutting can be included into the model 
(see also Painter [21]). This might necessitate to explicitly model the protease as a 
third variable through its own reaction-diffusion equation. 

A consideration of the length scales of the fibres relative to the size of the moving 
cells might also give valuable input into the appropriate modelling assumptions. 

Finally, we have studied an unbounded domain M" to avoid boundary conditions. 
To formulate the correct boundary conditions for model (7) is not trivial. A com- 
monly observed effect seen in tissue is that a tumour is encapsulated by a dense 
fibre network. In that case the fibres at the boundary will be aligned tangentially to 
the boundary and should trap moving cells inside the domain. The encapsulations 
can be understood as patchy steady states, as described here. A careful analysis of 
other boundary conditions and its implications on existence and steady states is left 
to future work. For the simulations in Figure 2 K. Painter used periodic boundary 
conditions on a square domain, i.e. a flat torus. 
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